  
# Load packages
  require(tidyverse)
  library(broom)
  library(lmtest)
  library(sandwich)
  library(stargazer)
  
# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]
  
# Prepare data
  load("ObsStudy_00_data_Sweden_2010.RData")
  dt <- x

  # Evaluation Political Parties
  colnames(dt)[seq(89, 95, 1)] <- c("rat_C_factor", "rat_M_factor", "rat_V_factor", "rat_FP_factor", "rat_S_factor", "rat_MP_factor", "rat_KD_factor")

  # Generate numeric party rating variables
  dt$rat_C <- ifelse(dt$rat_C_factor == "0 really dislike", 0, NA)
  dt$rat_C <- ifelse(dt$rat_C_factor == "1", 1, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "2", 2, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "3", 3, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "4", 4, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "5", 5, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "6", 6, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "7", 7, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "8", 8, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "9", 9, dt$rat_C)
  dt$rat_C <- ifelse(dt$rat_C_factor == "10 really like", 10, dt$rat_C)
  
  dt$rat_M <- ifelse(dt$rat_M_factor == "0 really dislike", 0, NA)
  dt$rat_M <- ifelse(dt$rat_M_factor == "1", 1, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "2", 2, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "3", 3, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "4", 4, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "5", 5, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "6", 6, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "7", 7, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "8", 8, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "9", 9, dt$rat_M)
  dt$rat_M <- ifelse(dt$rat_M_factor == "10 really like", 10, dt$rat_M)
  
  dt$rat_V <- ifelse(dt$rat_V_factor == "0 really dislike", 0, NA)
  dt$rat_V <- ifelse(dt$rat_V_factor == "1", 1, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "2", 2, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "3", 3, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "4", 4, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "5", 5, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "6", 6, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "7", 7, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "8", 8, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "9", 9, dt$rat_V)
  dt$rat_V <- ifelse(dt$rat_V_factor == "10 really like", 10, dt$rat_V)
  
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "0 really dislike", 0, NA)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "1", 1, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "2", 2, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "3", 3, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "4", 4, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "5", 5, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "6", 6, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "7", 7, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "8", 8, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "9", 9, dt$rat_FP)
  dt$rat_FP <- ifelse(dt$rat_FP_factor == "10 really like", 10, dt$rat_FP)
  
  dt$rat_S <- ifelse(dt$rat_S_factor == "0 really dislike", 0, NA)
  dt$rat_S <- ifelse(dt$rat_S_factor == "1", 1, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "2", 2, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "3", 3, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "4", 4, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "5", 5, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "6", 6, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "7", 7, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "8", 8, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "9", 9, dt$rat_S)
  dt$rat_S <- ifelse(dt$rat_S_factor == "10 really like", 10, dt$rat_S)
  
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "0 really dislike", 0, NA)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "1", 1, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "2", 2, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "3", 3, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "4", 4, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "5", 5, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "6", 6, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "7", 7, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "8", 8, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "9", 9, dt$rat_MP)
  dt$rat_MP <- ifelse(dt$rat_MP_factor == "10 really like", 10, dt$rat_MP)
  
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "0 really dislike", 0, NA)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "1", 1, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "2", 2, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "3", 3, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "4", 4, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "5", 5, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "6", 6, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "7", 7, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "8", 8, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "9", 9, dt$rat_KD)
  dt$rat_KD <- ifelse(dt$rat_KD_factor == "10 really like", 10, dt$rat_KD)
  
  
  # Coalition Evaluation
  colnames(dt)[seq(114, 117, 1)] <- c("rat_coal_M_FP_C_KD_factor", "rat_coal_S_MP_V_factor", "rat_coal_M_FP_C_KD_MP_factor", "rat_coal_S_MP_V_C_factor")

  # Generate numeric coalition evaluation variable
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "0 really dislike", 0, NA)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "1", 1, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "2", 2, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "3", 3, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "4", 4, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "5", 5, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "6", 6, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "7", 7, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "8", 8, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "9", 9, dt$rat_coal_M_FP_C_KD)
  dt$rat_coal_M_FP_C_KD <- ifelse(dt$rat_coal_M_FP_C_KD_factor == "10 really like", 10, dt$rat_coal_M_FP_C_KD)
  
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "0 really dislike", 0, NA)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "1", 1, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "2", 2, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "3", 3, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "4", 4, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "5", 5, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "6", 6, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "7", 7, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "8", 8, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "9", 9, dt$rat_coal_S_MP_V)
  dt$rat_coal_S_MP_V <- ifelse(dt$rat_coal_S_MP_V_factor == "10 really like", 10, dt$rat_coal_S_MP_V)
  
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "0 really dislike", 0, NA)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "1", 1, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "2", 2, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "3", 3, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "4", 4, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "5", 5, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "6", 6, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "7", 7, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "8", 8, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "9", 9, dt$rat_coal_M_FP_C_KD_MP)
  dt$rat_coal_M_FP_C_KD_MP <- ifelse(dt$rat_coal_M_FP_C_KD_MP_factor == "10 really like", 10, dt$rat_coal_M_FP_C_KD_MP)
  
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "0 really dislike", 0, NA)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "1", 1, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "2", 2, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "3", 3, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "4", 4, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "5", 5, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "6", 6, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "7", 7, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "8", 8, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "9", 9, dt$rat_coal_S_MP_V_C)
  dt$rat_coal_S_MP_V_C <- ifelse(dt$rat_coal_S_MP_V_C_factor == "10 really like", 10, dt$rat_coal_S_MP_V_C)
  
  
  # Coalition Likelihood
  colnames(dt)[seq(118, 121, 1)] <- c("coallik_M_FP_C_KD_factor", "coallik_S_MP_V_factor", "coallik_M_FP_C_KD_MP_factor", "coallik_S_MP_V_C_factor")

  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "0 not likely at all", 0, NA)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "1", 1, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "2", 2, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "3", 3, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "4", 4, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "5", 5, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "6", 6, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "7", 7, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "8", 8, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "9", 9, dt$coallik_M_FP_C_KD)
  dt$coallik_M_FP_C_KD <- ifelse(dt$coallik_M_FP_C_KD_factor == "10 very likely", 10, dt$coallik_M_FP_C_KD)
  
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "0 not likely at all", 0, NA)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "1", 1, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "2", 2, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "3", 3, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "4", 4, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "5", 5, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "6", 6, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "7", 7, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "8", 8, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "9", 9, dt$coallik_S_MP_V)
  dt$coallik_S_MP_V <- ifelse(dt$coallik_S_MP_V_factor == "10 very likely", 10, dt$coallik_S_MP_V)
  
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "0 not likely at all", 0, NA)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "1", 1, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "2", 2, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "3", 3, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "4", 4, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "5", 5, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "6", 6, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "7", 7, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "8", 8, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "9", 9, dt$coallik_M_FP_C_KD_MP)
  dt$coallik_M_FP_C_KD_MP <- ifelse(dt$coallik_M_FP_C_KD_MP_factor == "10 very likely", 10, dt$coallik_M_FP_C_KD_MP)
  
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "0 not likely at all", 0, NA)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "1", 1, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "2", 2, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "3", 3, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "4", 4, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "5", 5, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "6", 6, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "7", 7, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "8", 8, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "9", 9, dt$coallik_S_MP_V_C)
  dt$coallik_S_MP_V_C <- ifelse(dt$coallik_S_MP_V_C_factor == "10 very likely", 10, dt$coallik_S_MP_V_C)
  
  
  # Dependent Variable
  dt$vote <- ifelse(dt$pre1_Q15 == "left party", "V", NA)
  dt$vote <- ifelse(dt$pre1_Q15 == "social democrats", "S", dt$vote)
  dt$vote <- ifelse(dt$pre1_Q15 == "centre party", "C", dt$vote)
  dt$vote <- ifelse(dt$pre1_Q15 == "peoples liberal party", "FP", dt$vote)
  dt$vote <- ifelse(dt$pre1_Q15 == "moderate party", "M", dt$vote)
  dt$vote <- ifelse(dt$pre1_Q15 == "christian democrats", "KD", dt$vote)
  dt$vote <- ifelse(dt$pre1_Q15 == "green party", "MP", dt$vote)

  ## Generate ptv columns based on vote variable
  dt$ptv_V <- ifelse(dt$vote == "V", 1, 0)
  dt$ptv_S <- ifelse(dt$vote == "S", 1, 0)
  dt$ptv_C <- ifelse(dt$vote == "C", 1, 0)
  dt$ptv_FP <- ifelse(dt$vote == "FP", 1, 0)
  dt$ptv_M <- ifelse(dt$vote == "M", 1, 0)
  dt$ptv_KD <- ifelse(dt$vote == "KD", 1, 0)
  dt$ptv_MP <- ifelse(dt$vote == "MP", 1, 0)
  
  #### Coalition Evaluation
  Z <- dt %>% select("rat_coal_M_FP_C_KD", "rat_coal_S_MP_V", "rat_coal_M_FP_C_KD_MP", "rat_coal_S_MP_V_C") %>% mutate_all(funs(scales::rescale(.,to = c(0, 1), from = c(0,10)))) 
  
  rat_coal_M <- Z %>% select(rat_coal_M_FP_C_KD, rat_coal_M_FP_C_KD_MP) %>% as.matrix()
  rat_coal_FP <- Z %>% select(rat_coal_M_FP_C_KD, rat_coal_M_FP_C_KD_MP) %>% as.matrix()
  rat_coal_C <- Z %>% select(rat_coal_M_FP_C_KD, rat_coal_M_FP_C_KD_MP, rat_coal_S_MP_V_C) %>% as.matrix()
  rat_coal_KD <- Z %>% select(rat_coal_M_FP_C_KD, rat_coal_M_FP_C_KD_MP) %>% as.matrix()
  rat_coal_S <- Z %>% select(rat_coal_S_MP_V, rat_coal_S_MP_V_C) %>% as.matrix()
  rat_coal_MP <- Z %>% select(rat_coal_S_MP_V, rat_coal_M_FP_C_KD_MP, rat_coal_S_MP_V_C) %>% as.matrix()
  rat_coal_V <- Z %>% select(rat_coal_S_MP_V, rat_coal_S_MP_V_C) %>% as.matrix()
  
  
  ### Coalition Likelihood
  gamma <- dt %>% select("coallik_M_FP_C_KD", "coallik_S_MP_V", "coallik_M_FP_C_KD_MP", "coallik_S_MP_V_C") 
  
  coal_lik_M <- gamma %>% select(coallik_M_FP_C_KD, coallik_M_FP_C_KD_MP) %>%
    mutate(sum_exp = exp(coallik_M_FP_C_KD) + exp(coallik_M_FP_C_KD_MP),
           coallik_M_FP_C_KD = exp(coallik_M_FP_C_KD)/sum_exp,
           coallik_M_FP_C_KD_MP = exp(coallik_M_FP_C_KD_MP)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_FP <- gamma %>% select(coallik_M_FP_C_KD, coallik_M_FP_C_KD_MP) %>%
    mutate(sum_exp = exp(coallik_M_FP_C_KD) + exp(coallik_M_FP_C_KD_MP),
           coallik_M_FP_C_KD = exp(coallik_M_FP_C_KD)/sum_exp,
           coallik_M_FP_C_KD_MP = exp(coallik_M_FP_C_KD_MP)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_C <- gamma %>% select(coallik_M_FP_C_KD, coallik_M_FP_C_KD_MP, coallik_S_MP_V_C) %>%
    mutate(sum_exp = exp(coallik_M_FP_C_KD) + exp(coallik_M_FP_C_KD_MP) + exp(coallik_S_MP_V_C),
           coallik_M_FP_C_KD = exp(coallik_M_FP_C_KD)/sum_exp,
           coallik_M_FP_C_KD_MP = exp(coallik_M_FP_C_KD_MP)/sum_exp,
           coallik_S_MP_V_C = exp(coallik_S_MP_V_C)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_KD <- gamma %>% select(coallik_M_FP_C_KD, coallik_M_FP_C_KD_MP) %>%
    mutate(sum_exp = exp(coallik_M_FP_C_KD) + exp(coallik_M_FP_C_KD_MP),
           coallik_M_FP_C_KD = exp(coallik_M_FP_C_KD)/sum_exp,
           coallik_M_FP_C_KD_MP = exp(coallik_M_FP_C_KD_MP)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_S <- gamma %>% select(coallik_S_MP_V, coallik_S_MP_V_C) %>%
    mutate(sum_exp = exp(coallik_S_MP_V) + exp(coallik_S_MP_V_C),
           coallik_S_MP_V = exp(coallik_S_MP_V)/sum_exp,
           coallik_S_MP_V_C = exp(coallik_S_MP_V_C)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_MP <- gamma %>% select(coallik_S_MP_V, coallik_M_FP_C_KD_MP, coallik_S_MP_V_C) %>%
    mutate(sum_exp = exp(coallik_S_MP_V) + exp(coallik_M_FP_C_KD_MP) + exp(coallik_S_MP_V_C),
           coallik_S_MP_V = exp(coallik_S_MP_V)/sum_exp,
           coallik_M_FP_C_KD_MP = exp(coallik_M_FP_C_KD_MP)/sum_exp,
           coallik_S_MP_V_C = exp(coallik_S_MP_V_C)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_V <- gamma %>% select(coallik_S_MP_V, coallik_S_MP_V_C) %>%
    mutate(sum_exp = exp(coallik_S_MP_V) + exp(coallik_S_MP_V_C),
           coallik_S_MP_V = exp(coallik_S_MP_V)/sum_exp,
           coallik_S_MP_V_C = exp(coallik_S_MP_V_C)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  # Mean-Variance Model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  M_M <- sapply(1:N, function(i) EV(coal_lik_M[i,], rat_coal_M[i,]))
  V_M <- sapply(1:N, function(i) Vari(coal_lik_M[i,], rat_coal_M[i,]))
  
  M_FP <- sapply(1:N, function(i) EV(coal_lik_FP[i,], rat_coal_FP[i,]))
  V_FP <- sapply(1:N, function(i) Vari(coal_lik_FP[i,], rat_coal_FP[i,]))
  
  M_C <- sapply(1:N, function(i) EV(coal_lik_C[i,], rat_coal_C[i,]))
  V_C <- sapply(1:N, function(i) Vari(coal_lik_C[i,], rat_coal_C[i,]))
  
  M_KD <- sapply(1:N, function(i) EV(coal_lik_KD[i,], rat_coal_KD[i,]))
  V_KD <- sapply(1:N, function(i) Vari(coal_lik_KD[i,], rat_coal_KD[i,]))
  
  M_S <- sapply(1:N, function(i) EV(coal_lik_S[i,], rat_coal_S[i,]))
  V_S <- sapply(1:N, function(i) Vari(coal_lik_S[i,], rat_coal_S[i,]))
  
  M_MP <- sapply(1:N, function(i) EV(coal_lik_MP[i,], rat_coal_MP[i,]))
  V_MP <- sapply(1:N, function(i) Vari(coal_lik_MP[i,], rat_coal_MP[i,]))
  
  M_V <- sapply(1:N, function(i) EV(coal_lik_V[i,], rat_coal_V[i,]))
  V_V <- sapply(1:N, function(i) Vari(coal_lik_V[i,], rat_coal_V[i,]))
  
  # Gender
  dt$male <- ifelse(dt$r_sex == "man", 1, 0)

  # Age
  dt$age <- ifelse(dt$r_age >=18, dt$r_age, NA)

  # Education
  dt$edu <- ifelse(dt$post2_Q14 == "not completed primary or comprehensive school", 1, NA)
  dt$edu <- ifelse(dt$post2_Q14 == "primary or comprehensive school", 2, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "studies at highschool or equivalent", 3, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "graduation from highschool or equivalent", 4, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "post-high school studies, not university", 5, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "studies at university", 6, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "graduation from university", 7, dt$edu)
  dt$edu <- ifelse(dt$post2_Q14 == "graduation from post-graduate studies", 8, dt$edu)
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$pid_M <- ifelse(dt$pre1_Q31=="the Moderate party", 1, 0)
  dt$pid_FP <- ifelse(dt$pre1_Q31=="the Peoples liberal party", 1, 0)
  dt$pid_C <- ifelse(dt$pre1_Q31=="the Centre party", 1, 0)
  dt$pid_KD <- ifelse(dt$pre1_Q31=="the Christian democrats", 1, 0)
  dt$pid_S <- ifelse(dt$pre1_Q31=="the Social democrats", 1, 0)
  dt$pid_MP <- ifelse(dt$pre1_Q31=="the Green party", 1, 0)
  dt$pid_V <- ifelse(dt$pre1_Q31=="the Left party", 1, 0)

  d <- data.frame(
    "ptv_M" = dt$ptv_M,
    "ptv_FP" = dt$ptv_FP,
    "ptv_C" = dt$ptv_C,
    "ptv_KD" = dt$ptv_KD,
    "ptv_S" = dt$ptv_S,
    "ptv_MP" = dt$ptv_MP,
    "ptv_V" = dt$ptv_V,
    "rat.M" = dt$rat_M ,
    "rat.FP" = dt$rat_FP,
    "rat.C" = dt$rat_C,
    "rat.KD" = dt$rat_KD,
    "rat.S" = dt$rat_S,
    "rat.MP" = dt$rat_MP,
    "rat.V" = dt$rat_V,
    "pid_M" = dt$pid_M,
    "pid_FP" = dt$pid_FP,
    "pid_C" = dt$pid_C,
    "pid_KD" = dt$pid_KD,
    "pid_S" = dt$pid_S,
    "pid_MP" = dt$pid_MP,
    "pid_V" = dt$pid_V,
    "lotterymean.M" = M_M,
    "lotteryvariance.M" = V_M,
    "lotterymean.FP" = M_FP,
    "lotteryvariance.FP" = V_FP,
    "lotterymean.C" = M_C,
    "lotteryvariance.C" = V_C,
    "lotterymean.KD" = M_KD,
    "lotteryvariance.KD" = V_KD,
    "lotterymean.S" = M_S,
    "lotteryvariance.S" = V_S,
    "lotterymean.MP" = M_MP,
    "lotteryvariance.MP" = V_MP,
    "lotterymean.V" = M_V,
    "lotteryvariance.V" = V_V,
    "sex" = dt$male,
    "age" = dt$age,
    "edu" = dt$edu
  )
  
  d$lotteryvariance.M <- scales::rescale(d$lotteryvariance.M, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.FP <- scales::rescale(d$lotteryvariance.FP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.C <- scales::rescale(d$lotteryvariance.C, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.KD <- scales::rescale(d$lotteryvariance.KD, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.S <- scales::rescale(d$lotteryvariance.S, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.MP <- scales::rescale(d$lotteryvariance.MP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.V <- scales::rescale(d$lotteryvariance.V, to = c(0, 1), from = c(0,0.25))
  
# Saving model results
  
  m1_tidy <- tidy(m1 <- lm(as.formula(paste("ptv_M ~ lotteryvariance.M + lotterymean.M + rat.M + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_M" else "")), d))
  m2_tidy <- tidy(m2 <- lm(as.formula(paste("ptv_FP ~ lotteryvariance.FP + lotterymean.FP + rat.FP + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_FP" else "")), d))
  m3_tidy <- tidy(m3 <- lm(as.formula(paste("ptv_C ~ lotteryvariance.C + lotterymean.C + rat.C + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_C" else "")), d))
  m4_tidy <- tidy(m4 <- lm(as.formula(paste("ptv_KD ~ lotteryvariance.KD + lotterymean.KD + rat.KD + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_KD" else "")), d))
  m5_tidy <- tidy(m5 <- lm(as.formula(paste("ptv_S ~ lotteryvariance.S + lotterymean.S + rat.S + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_S" else "")), d))
  m6_tidy <- tidy(m6 <- lm(as.formula(paste("ptv_MP ~ lotteryvariance.MP + lotterymean.MP + rat.MP + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_MP" else "")), d))
  m7_tidy <- tidy(m7 <- lm(as.formula(paste("ptv_V ~ lotteryvariance.V + lotterymean.V + rat.V + sex + as.factor(edu) + age", 
                                            if (robustnesscheck_pid) "+ pid_V" else "")), d))
  
  SE_2010_ICP_M_Variance_Estimate <- m1_tidy %>% filter(term=="lotteryvariance.M") %>% select("estimate")
  SE_2010_ICP_M_Variance_SE <- se_robust(m1)["lotteryvariance.M"] #m1_tidy %>% filter(term=="lotteryvariance.M") %>% select("std.error")
  SE_2010_ICP_M_Mean_Estimate <- m1_tidy %>% filter(term=="lotterymean.M") %>% select("estimate")
  SE_2010_ICP_M_Mean_SE <- se_robust(m1)["lotterymean.M"] #m1_tidy %>% filter(term=="lotterymean.M") %>% select("std.error")
  
  SE_2010_ICP_FP_Variance_Estimate <- m2_tidy %>% filter(term=="lotteryvariance.FP") %>% select("estimate")
  SE_2010_ICP_FP_Variance_SE <- se_robust(m2)["lotteryvariance.FP"] #m2_tidy %>% filter(term=="lotteryvariance.FP") %>% select("std.error")
  SE_2010_ICP_FP_Mean_Estimate <- m2_tidy %>% filter(term=="lotterymean.FP") %>% select("estimate")
  SE_2010_ICP_FP_Mean_SE <- se_robust(m2)["lotterymean.FP"] #m2_tidy %>% filter(term=="lotterymean.FP") %>% select("std.error")
  
  SE_2010_ICP_C_Variance_Estimate <- m3_tidy %>% filter(term=="lotteryvariance.C") %>% select("estimate")
  SE_2010_ICP_C_Variance_SE <- se_robust(m3)["lotteryvariance.C"] #m3_tidy %>% filter(term=="lotteryvariance.C") %>% select("std.error")
  SE_2010_ICP_C_Mean_Estimate <- m3_tidy %>% filter(term=="lotterymean.C") %>% select("estimate")
  SE_2010_ICP_C_Mean_SE <- se_robust(m3)["lotterymean.C"] #m3_tidy %>% filter(term=="lotterymean.C") %>% select("std.error")
  
  SE_2010_ICP_KD_Variance_Estimate <- m4_tidy %>% filter(term=="lotteryvariance.KD") %>% select("estimate")
  SE_2010_ICP_KD_Variance_SE <- se_robust(m4)["lotteryvariance.KD"] #m4_tidy %>% filter(term=="lotteryvariance.KD") %>% select("std.error")
  SE_2010_ICP_KD_Mean_Estimate <- m4_tidy %>% filter(term=="lotterymean.KD") %>% select("estimate")
  SE_2010_ICP_KD_Mean_SE <- se_robust(m4)["lotterymean.KD"] #m4_tidy %>% filter(term=="lotterymean.KD") %>% select("std.error")
  
  SE_2010_ICP_S_Variance_Estimate <- m5_tidy %>% filter(term=="lotteryvariance.S") %>% select("estimate")
  SE_2010_ICP_S_Variance_SE <- se_robust(m5)["lotteryvariance.S"] #m5_tidy %>% filter(term=="lotteryvariance.S") %>% select("std.error")
  SE_2010_ICP_S_Mean_Estimate <- m5_tidy %>% filter(term=="lotterymean.S") %>% select("estimate")
  SE_2010_ICP_S_Mean_SE <- se_robust(m5)["lotterymean.S"] #m5_tidy %>% filter(term=="lotterymean.S") %>% select("std.error")
  
  SE_2010_ICP_MP_Variance_Estimate <- m6_tidy %>% filter(term=="lotteryvariance.MP") %>% select("estimate")
  SE_2010_ICP_MP_Variance_SE <- se_robust(m6)["lotteryvariance.MP"] #m6_tidy %>% filter(term=="lotteryvariance.MP") %>% select("std.error")
  SE_2010_ICP_MP_Mean_Estimate <- m6_tidy %>% filter(term=="lotterymean.MP") %>% select("estimate")
  SE_2010_ICP_MP_Mean_SE <- se_robust(m6)["lotterymean.MP"] #m6_tidy %>% filter(term=="lotterymean.MP") %>% select("std.error")
  
  SE_2010_ICP_V_Variance_Estimate <- m7_tidy %>% filter(term=="lotteryvariance.V") %>% select("estimate")
  SE_2010_ICP_V_Variance_SE <- se_robust(m7)["lotteryvariance.V"] #m7_tidy %>% filter(term=="lotteryvariance.V") %>% select("std.error")
  SE_2010_ICP_V_Mean_Estimate <- m7_tidy %>% filter(term=="lotterymean.V") %>% select("estimate")
  SE_2010_ICP_V_Mean_SE <- se_robust(m7)["lotterymean.V"] #m7_tidy %>% filter(term=="lotterymean.V") %>% select("std.error")

  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.M"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.M"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.FP"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.FP"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.C"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.C"] <- "Government Lottery Mean"
  names(m4$coefficients)[names(m4$coefficients) == "lotteryvariance.KD"] <- "Government Lottery Variance"
  names(m4$coefficients)[names(m4$coefficients) == "lotterymean.KD"] <- "Government Lottery Mean"
  names(m5$coefficients)[names(m5$coefficients) == "lotteryvariance.S"] <- "Government Lottery Variance"
  names(m5$coefficients)[names(m5$coefficients) == "lotterymean.S"] <- "Government Lottery Mean"
  names(m6$coefficients)[names(m6$coefficients) == "lotteryvariance.MP"] <- "Government Lottery Variance"
  names(m6$coefficients)[names(m6$coefficients) == "lotterymean.MP"] <- "Government Lottery Mean"
  names(m7$coefficients)[names(m7$coefficients) == "lotteryvariance.V"] <- "Government Lottery Variance"
  names(m7$coefficients)[names(m7$coefficients) == "lotterymean.V"] <- "Government Lottery Mean"
  
  if(!robustnesscheck_pid){
  
  stargazer(m1, m2, m3, m4, m5, m6, m7, title="Linear regressions of vote choice on perceived government lottery variance and mean (The 2010 Internet Campaign Panel, Sweden)", align=TRUE, 
            dep.var.labels=c("M", "FP", "C", "KD", "S", "MP", "V"),
            omit.stat=c("LL","ser","f"),
            se = lapply(list(m1, m2, m3, m4, m5, m6, m7), se_robust),
            p = lapply(list(m1, m2, m3, m4, m5, m6, m7), p_robust),
            type = "text",
            out = paste0("TableSM15.tex")
            )

  }